package body ariane.maths.matrix.square is

  -- creates an object and initializes it with specified data
  function create(components : components_t := (others=>(others=>zero)))
                  return object is
    o : object;
  begin
    o.set(components);
    return o;
  end create;

  -- constructor of matrix object
  -- as matrix type here is a constrained, we can use procedure-based
  -- initializer-styled constructor for it
  -- the component array argument is optional which when absent is
  -- an array of default values
  procedure set(o : in out object;
                components : components_t := (others=>(others=>zero))) is
  begin
    o.components := components;
  end set;

  -- constructor of vector object
  -- the object is initialized as an identity matrix
  procedure identify(o : in out object) is
  begin
    o.components := (others=>(others=>zero));
    for i in 1..dimension loop
      o.components(i,i) := one;
    end loop;
  end identify;

  -- setter that assigns value to specific component
  procedure set(o : in out object; row, col : positive; val : component_t) is
  begin
    o.components(row, col) := val;
  end set;

  -- getter that returns value of specific component
  function get(o : in object; row, col : positive) return component_t is
  begin
    return o.components(row, col);
  end get;

  -- returns the result of adding matrix lhs to rhs
  function "+"(lhs, rhs: object) return object is
    result : object;
  begin
    for i in 1..dimension loop
      for j in 1..dimension loop
        result.components(i,j) := lhs.components(i,j) + rhs.components(i,j);
      end loop;
    end loop;
    return result;
  end "+";

  -- returns the result of subtracting matrix rhs from lhs
  function "-"(lhs, rhs: object) return object is
    result : object;
  begin
    for i in 1..dimension loop
      for j in 1..dimension loop
        result.components(i,j) := lhs.components(i,j) - rhs.components(i,j);
      end loop;
    end loop;
    return result;
  end "-";

  -- returns the product of two matrices
  -- todo : optimization might be possible
  function "*"(lhs, rhs: object) return object is
    result : object;
  begin
    for i in 1..dimension loop
      for j in 1..dimension loop
        result.components(i,j) := zero;
        for k in 1..dimension loop
          result.components(i,j) := result.components(i,j) +
            lhs.components(i,k) * rhs.components(k,j);
        end loop;
      end loop;
    end loop;
    return result;
  end "*";

  -- returns the negative of the given object
  function "-"(o : object) return object is
    result : object;
  begin
    for i in 1..dimension loop
      for j in 1..dimension loop
        result.components(i,j) := -o.components(i,j);
      end loop;
    end loop;
    return result;
  end "-";

  -- returns the result of a matrix right multiplied by a vector
  function "*"(o : object; v : vector.object) return vector.object is
    tmp : component_t;
    res : vector.object;
  begin
    for i in 1..dimension loop
      tmp := zero;
      for j in 1..dimension loop
        tmp := tmp + o.components(i,j) * v.get(j);
      end loop;
      res.set(i, tmp);
    end loop;
    return res;
  end "*";

  -- returns the result of a matrix left multiplied by a vector
  function "*"(v : vector.object; o : object) return vector.object is
    tmp : component_t;
    res : vector.object;
  begin
    for i in 1..dimension loop
      tmp := zero;
      for j in 1..dimension loop
        tmp := tmp + v.get(j) * o.components(j,i);
      end loop;
      res.set(i, tmp);
    end loop;
    return res;
  end "*";

  -- turns the matrix to its inverse
  -- this uses guassian method
  procedure invert(o : in out object; inverted : out boolean) is

    partner : object;	-- matrix that acompanies this in transformation

    -- swap two rows indexed by a and b in both matrices
    procedure swap(a, b : positive) is
      tmp : component_t;
    begin
      for i in 1..dimension loop
        tmp := o.components(a,i);
        o.components(a,i) := o.components(b,i);
        o.components(b,i) := tmp;
        tmp := partner.components(a,i);
        partner.components(a,i) := partner.components(b,i);
        partner.components(b,i) := tmp;
      end loop;
    end swap;

    -- exchange the dominant row (row with largest leading component
    -- starting i stage the row to the first row to process against
    -- true return if a dominant row with nonzero leading component is found
    function pickdominantrow(i : positive) return boolean is
      -- variables for locating row with largest leading component
      -- largest leading component
      largestleading : component_t:= abs o.components(i,i);
      rowll : positive := i; -- row with largest leading element
    begin
      -- pick the row with largest leading component to process against
      for j in i+1..dimension loop
        if largestleading < abs o.components(j,i) then
          largestleading := abs o.components(j,i);
          rowll := j;
        end if;
      end loop;

      if largestleading = zero then
        return false;
      end if;

      if rowll /= i then
        swap(i,rowll);
      end if;

      return true;
    end pickdominantrow;

    -- normalize the dominant row by dividing it by its leading component
    -- idxll is the index of the lifted dominant row as well as
    -- of the leading component
    procedure normalizedominant(idxll : positive) is
      -- largest leading component to divide the row by
      ll : component_t := o.components(idxll,idxll);
    begin
      -- todo check and see if it is worth optimizing by turning division
      -- to multiplication
      for i in idxll+1..dimension loop
        o.components(idxll,i) := o.components(idxll,i) / ll;
      end loop;
      -- apply the row transformation to the accompanying matrix
      for i in 1..dimension loop
        partner.components(idxll,i) := partner.components(idxll,i) / ll;
      end loop;
      o.components(idxll,idxll) := one;
    end normalizedominant;

    -- zero the leading component of subordinate rows according to dominant
    -- idxll is the index of dominant
    -- dominant row has already been normalized
    procedure beheadsubordinate(idxll : positive) is
      lc : component_t; -- leading component of each row
    begin
      for i in idxll+1..dimension loop
        if o.components(i,idxll) /= zero then
          lc := o.components(i,idxll);
          o.components(i,idxll) := zero;
          for j in idxll+1..dimension loop
            o.components(i,j) := o.components(i,j) - o.components(idxll,j) * lc;
          end loop;
          for j in 1..dimension loop
            partner.components(i,j) := partner.components(i,j) -
              partner.components(idxll,j) * lc;
          end loop;
        end if;
      end loop;
    end beheadsubordinate;

    -- remove everything except the leading one from the normalized row
    -- same row transformation is applied to the accompanying matrix
    procedure eviscerate(row : positive) is
    begin
      for j in 1..dimension loop
        for k in row+1..dimension loop
          partner.components(row,j) := partner.components(row,j)
            - partner.components(k,j) * o.components(row,k);
        end loop;
      end loop;
      for j in row+1..dimension loop
        o.components(row,j) := zero;
      end loop;
    end eviscerate;

  begin
    partner.identify;	-- partner is identity initally

    for i in 1..dimension loop
      if not pickdominantrow(i) then
        inverted := false;
        return;
      end if;
      normalizedominant(i);
      beheadsubordinate(i);
    end loop;

    for i in reverse 1..dimension-1 loop
      eviscerate(i);
    end loop;

    inverted := true;
    o := partner;	-- partner is the inverse of the original matrix
  end invert;

  -- returns the inverse of the give matrix
  -- todo implement it using some approach
  function inverse(o : in object) return object is
    result : object := o;
    inverted : boolean;
  begin
    result.invert(inverted);
    if not inverted then
      result.set; -- return a zero matrix to indicate non-invertible
    end if;
    return result;
  end inverse;

    -- gets the number of rows in the matrix
  function numrows(o : in object) return positive is
  begin
    return dimension;
  end numrows;

  -- gets the number of colums in the matrix
  function numcols(o : in object) return positive is
  begin
    return dimension;
  end numcols;


end ariane.maths.matrix.square;
